{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "using VMLS\n",
    "using LinearAlgebra"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Chapter 10\n",
    "# Matrix multiplication\n",
    "### 10.1 Matrix-matrix multiplication\n",
    "In Julia the product of matrices `A` and `B` is obtained with `A*B`. We calculate the matrix product on page [177](https://web.stanford.edu/~boyd/vmls/vmls.pdf#section.10.1) of VMLS. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2×3 Array{Float64,2}:\n",
       " -1.5   3.0  2.0\n",
       "  1.0  -1.0  0.0"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A = [-1.5 3 2; 1 -1 0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3×2 Array{Int64,2}:\n",
       " -1  -1\n",
       "  0  -2\n",
       "  1   0"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "B = [-1 -1; 0 -2; 1 0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2×2 Array{Float64,2}:\n",
       "  3.5  -4.5\n",
       " -1.0   1.0"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "C = A*B"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Gram matrix.** The Gram matrix of a matrix $A$ is the matrix $G = A^TA$. It is a symmetric matrix and the $i$, $j$ element $G_{ij}$ is the inner product of columns $i$ and $j$ of $A$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3×3 Array{Float64,2}:\n",
       "  9.09318  -3.07972   2.43275\n",
       " -3.07972  11.1709   -4.92137\n",
       "  2.43275  -4.92137  13.0071 "
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    " A = randn(10,3);\n",
    " G = A'*A"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "11.1709135648355"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Gii is norm of column i, squared\n",
    "G[2,2]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "11.170913564835498"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "norm(A[:,2])^2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(2.432746740923237, 2.4327467409232364)"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Gij is inner product of columns i and j\n",
    "G[1,3],A[:,1]'*A[:,3]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Complexity of matrix triple product.** Let’s check the associative property, which\n",
    "states that $(AB)C = A(BC)$ for any $m×n$ matrix $A$, any $n×p$ matrix $B$, and any\n",
    "$p × q$ matrix $B$. At the same time we will see that the left-hand and right-hand\n",
    "sides take very different amounts of time to compute."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "  0.124467 seconds (8 allocations: 61.035 MiB, 0.84% gc time)\n",
      "  0.169644 seconds (8 allocations: 61.035 MiB, 32.29% gc time)\n",
      "  0.018569 seconds (8 allocations: 31.281 MiB, 35.17% gc time)\n",
      "  0.013540 seconds (8 allocations: 31.281 MiB, 21.50% gc time)\n",
      "  0.126136 seconds (8 allocations: 61.035 MiB, 2.68% gc time)\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "5.298113387113036e-10"
      ]
     },
     "execution_count": 43,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "m = 2000; n = 50; q = 2000; p = 2000;\n",
    "A = randn(m,n); B = randn(n,p); C = randn(p,q);\n",
    "@time (A*B)*C\n",
    "@time (A*B)*C\n",
    "@time A*(B*C)\n",
    "@time A*(B*C)\n",
    "@time D = A*B*C; # evaluated as (A*B)*C or as A*(B*C)?\n",
    "\n",
    "LHS = (A*B)*C\n",
    "RHS = A*(B*C)\n",
    "norm(LHS-RHS)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We see that evaluating `(A*B)*C` takes around 10 times as much time as evaluating `A*(B*C)`, which is predicted from the complexities. In the last line we deduce that `A*B*C` is evaluated left to right, as `(A*B)*C`. Note that for these particular matrices, this is the (much) slower order to multiply the matrices.\n",
    "\n",
    "### 10.2 Composition of linear functions\n",
    "**Second difference matrix.** We compute the second difference matrix on page [184](https://web.stanford.edu/\\\\%7Eboyd/vmls/vmls.pdf#section*.221) of VMLS.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 51,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "4×5 Array{Float64,2}:\n",
       " -1.0   1.0   0.0   0.0  0.0\n",
       "  0.0  -1.0   1.0   0.0  0.0\n",
       "  0.0   0.0  -1.0   1.0  0.0\n",
       "  0.0   0.0   0.0  -1.0  1.0"
      ]
     },
     "execution_count": 51,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "d(n) = [-eye(n-1) zeros(n-1)] + [zeros(n-1) eye(n-1)];\n",
    "d(5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 52,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3×5 Array{Float64,2}:\n",
       " 1.0  -2.0   1.0   0.0  0.0\n",
       " 0.0   1.0  -2.0   1.0  0.0\n",
       " 0.0   0.0   1.0  -2.0  1.0"
      ]
     },
     "execution_count": 52,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "d(4)\n",
    "Delta = d(4)*d(5) # Second difference matrix"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 10.3 Matrix power\n",
    "The $k$th power of a square matrix $A$ is denoted $A^k$. In Julia, this power is formed\n",
    "using `A^k`.\n",
    "\n",
    "Let’s form the adjacency matrix of the directed graph on VMLS page [186](https://web.stanford.edu/~boyd/vmls/vmls.pdf#section*.226). Then let’s find out how many cycles of length $8$ there are, starting from each node. (A cycle is a path that starts and stops at the same node.) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 54,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "5×5 Array{Int64,2}:\n",
       " 0  1  0  0  1\n",
       " 1  0  1  0  0\n",
       " 0  0  1  1  1\n",
       " 1  0  0  0  0\n",
       " 0  0  0  1  0"
      ]
     },
     "execution_count": 54,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A = [ 0 1 0 0 1; 1 0 1 0 0; 0 0 1 1 1; 1 0 0 0 0; 0 0 0 1 0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 55,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "5×5 Array{Int64,2}:\n",
       " 1  0  1  1  0\n",
       " 0  1  1  1  2\n",
       " 1  0  1  2  1\n",
       " 0  1  0  0  1\n",
       " 1  0  0  0  0"
      ]
     },
     "execution_count": 55,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A^2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 56,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "5×5 Array{Int64,2}:\n",
       " 18  11  15  20  20\n",
       " 25  14  21  28  26\n",
       " 24  14  20  27  26\n",
       " 11   6   9  12  11\n",
       "  6   4   5   7   7"
      ]
     },
     "execution_count": 56,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A^8"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 57,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "5-element Array{Int64,1}:\n",
       " 18\n",
       " 14\n",
       " 20\n",
       " 12\n",
       "  7"
      ]
     },
     "execution_count": 57,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "number_of_cycles = diag(A^8)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Population dynamics.** Let’s write the code for figure [10.2](https://web.stanford.edu/\\\\%7Eboyd/vmls/vmls.pdf#figure.10.2) in VMLS, which plots the contribution factor to the total US population in $2020$ (ignoring immigration), for each age in $2010$. The Julia plot is in figure 10.1. We can see that, not surprisingly, $20–25$ year olds have the highest contributing factor, around 1.5. This means that on average, each $20-25$ year old in $2010$ will be responsible for around $1.5$ people in $2020$. This takes into account any children they may have before then, and (as a very small effect) the few of them who will no longer be with us in $2020$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 59,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/svg+xml": [
       "<?xml version=\"1.0\" encoding=\"utf-8\"?>\n",
       "<svg xmlns=\"http://www.w3.org/2000/svg\" xmlns:xlink=\"http://www.w3.org/1999/xlink\" width=\"600\" height=\"400\" viewBox=\"0 0 2400 1600\">\n",
       "<defs>\n",
       "  <clipPath id=\"clip0200\">\n",
       "    <rect x=\"0\" y=\"0\" width=\"2400\" height=\"1600\"/>\n",
       "  </clipPath>\n",
       "</defs>\n",
       "<path clip-path=\"url(#clip0200)\" d=\"\n",
       "M0 1600 L2400 1600 L2400 0 L0 0  Z\n",
       "  \" fill=\"#ffffff\" fill-rule=\"evenodd\" fill-opacity=\"1\"/>\n",
       "<defs>\n",
       "  <clipPath id=\"clip0201\">\n",
       "    <rect x=\"480\" y=\"0\" width=\"1681\" height=\"1600\"/>\n",
       "  </clipPath>\n",
       "</defs>\n",
       "<path clip-path=\"url(#clip0200)\" d=\"\n",
       "M215.754 1425.62 L2352.76 1425.62 L2352.76 47.2441 L215.754 47.2441  Z\n",
       "  \" fill=\"#ffffff\" fill-rule=\"evenodd\" fill-opacity=\"1\"/>\n",
       "<defs>\n",
       "  <clipPath id=\"clip0202\">\n",
       "    <rect x=\"215\" y=\"47\" width=\"2138\" height=\"1379\"/>\n",
       "  </clipPath>\n",
       "</defs>\n",
       "<polyline clip-path=\"url(#clip0202)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  255.871,1425.62 255.871,47.2441 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip0202)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  764.972,1425.62 764.972,47.2441 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip0202)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  1274.07,1425.62 1274.07,47.2441 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip0202)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  1783.17,1425.62 1783.17,47.2441 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip0202)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  2292.27,1425.62 2292.27,47.2441 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip0202)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  215.754,1386.61 2352.76,1386.61 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip0202)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  215.754,952.066 2352.76,952.066 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip0202)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  215.754,517.524 2352.76,517.524 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip0202)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  215.754,82.9819 2352.76,82.9819 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip0200)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  215.754,1425.62 2352.76,1425.62 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip0200)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  215.754,1425.62 215.754,47.2441 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip0200)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  255.871,1425.62 255.871,1404.94 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip0200)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  764.972,1425.62 764.972,1404.94 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip0200)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  1274.07,1425.62 1274.07,1404.94 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip0200)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  1783.17,1425.62 1783.17,1404.94 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip0200)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  2292.27,1425.62 2292.27,1404.94 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip0200)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  215.754,1386.61 247.809,1386.61 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip0200)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  215.754,952.066 247.809,952.066 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip0200)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  215.754,517.524 247.809,517.524 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip0200)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  215.754,82.9819 247.809,82.9819 \n",
       "  \"/>\n",
       "<g clip-path=\"url(#clip0200)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:middle;\" transform=\"rotate(0, 255.871, 1479.62)\" x=\"255.871\" y=\"1479.62\">0</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip0200)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:middle;\" transform=\"rotate(0, 764.972, 1479.62)\" x=\"764.972\" y=\"1479.62\">25</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip0200)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:middle;\" transform=\"rotate(0, 1274.07, 1479.62)\" x=\"1274.07\" y=\"1479.62\">50</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip0200)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:middle;\" transform=\"rotate(0, 1783.17, 1479.62)\" x=\"1783.17\" y=\"1479.62\">75</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip0200)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:middle;\" transform=\"rotate(0, 2292.27, 1479.62)\" x=\"2292.27\" y=\"1479.62\">100</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip0200)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:end;\" transform=\"rotate(0, 191.754, 1404.11)\" x=\"191.754\" y=\"1404.11\">0.0</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip0200)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:end;\" transform=\"rotate(0, 191.754, 969.566)\" x=\"191.754\" y=\"969.566\">0.5</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip0200)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:end;\" transform=\"rotate(0, 191.754, 535.024)\" x=\"191.754\" y=\"535.024\">1.0</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip0200)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:end;\" transform=\"rotate(0, 191.754, 100.482)\" x=\"191.754\" y=\"100.482\">1.5</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip0200)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:66px; text-anchor:middle;\" transform=\"rotate(0, 1284.25, 1559.48)\" x=\"1284.25\" y=\"1559.48\">Age</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip0200)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:66px; text-anchor:middle;\" transform=\"rotate(-90, 89.2861, 736.431)\" x=\"89.2861\" y=\"736.431\">Factor</text>\n",
       "</g>\n",
       "<polyline clip-path=\"url(#clip0202)\" style=\"stroke:#009af9; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  276.235,524.354 296.599,518.861 316.963,518.42 337.327,518.126 357.691,517.91 378.055,517.772 398.419,503.042 418.784,488.517 439.148,474.087 459.512,459.802 \n",
       "  479.876,445.597 500.24,407.452 520.604,369.546 540.968,331.675 561.332,293.823 581.696,255.994 602.06,224.79 622.424,193.545 642.788,162.18 663.152,130.665 \n",
       "  683.516,99.0764 703.88,96.5947 724.244,94.0006 744.608,91.4204 764.972,88.8682 785.336,86.2547 805.7,113.783 826.064,141.26 846.428,168.836 866.792,196.444 \n",
       "  887.156,224.146 907.52,261.769 927.884,299.469 948.248,337.329 968.612,375.316 988.976,413.564 1009.34,433.77 1029.7,454.209 1050.07,474.843 1070.43,495.639 \n",
       "  1090.8,516.677 1111.16,523.176 1131.52,529.825 1151.89,536.676 1172.25,543.729 1192.62,550.876 1212.98,554.052 1233.34,557.453 1253.71,560.952 1274.07,564.713 \n",
       "  1294.44,568.653 1314.8,572.324 1335.17,576.234 1355.53,580.339 1375.89,584.943 1396.26,590.417 1416.62,595.864 1436.99,601.707 1457.35,608.071 1477.71,615.273 \n",
       "  1498.08,623.067 1518.44,631.592 1538.81,640.745 1559.17,650.909 1579.53,661.777 1599.9,673.154 1620.26,686.17 1640.63,700.641 1660.99,715.937 1681.35,732.575 \n",
       "  1701.72,750.017 1722.08,769.478 1742.45,790.956 1762.81,814.147 1783.17,839.063 1803.54,866.514 1823.9,895.911 1844.27,926.905 1864.63,959.966 1884.99,994.125 \n",
       "  1905.36,1030.83 1925.72,1069.34 1946.09,1106.12 1966.45,1143.75 1986.81,1180.11 2007.18,1213.69 2027.54,1245.15 2047.91,1273.11 2068.27,1296.74 2088.63,1317.89 \n",
       "  2109,1386.61 2129.36,1386.61 2149.73,1386.61 2170.09,1386.61 2190.45,1386.61 2210.82,1386.61 2231.18,1386.61 2251.55,1386.61 2271.91,1386.61 2292.27,1386.61 \n",
       "  \n",
       "  \"/>\n",
       "</svg>\n"
      ]
     },
     "execution_count": 59,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "D = population_data();\n",
    "birth = D[\"birth_rate\"];\n",
    "death = D[\"death_rate\"];\n",
    "# Dynamics matrix for populaion dynamics\n",
    "A = [birth'; diagonal(1 .- death[1:end-1]) zeros(length(death)-1)];\n",
    "# Contribution factor to total poulation in 2020\n",
    "# from each age in 2010\n",
    "cf = ones(100)'*(A^10); # Contribution factor\n",
    "using Plots\n",
    "plot(cf', legend = false, xlabel = \"Age\", ylabel = \"Factor\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Figure 10.1** Contribution factor per age in $2010$ to the total population in $2020$. The value for age i− 1 is the ith component of the row vector 1TA10.\n",
    "\n",
    "### 10.4 QR factorization\n",
    "\n",
    "In Julia, the $QR$ factorization of a matrix $A$ can be found using `qr(A)`, which\n",
    "returns a tuple with the $Q$ and $R$ factors. However the matrix $Q$ is not returned as an array, but in a special compact format. It can be converted to a regular matrix\n",
    "variable using the command `Matrix(Q)`. Hence, the `QR` factorization as defined in\n",
    "VMLS is computed by a sequence of two commands:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "```julia\n",
    "Q,R = qr(A)\n",
    "Q = Matrix(Q)\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following example also illustates a second, but minor difference with the\n",
    "VMLS definition. The $R$ factor computed by Julia may have negative elements on\n",
    "the diagonal, as opposed to only positive elements if we follow the definition used\n",
    "in VMLS. The two definitions are equivalent, because if $R_{ii}$ is negative, one can\n",
    "change the sign of the ith row of R and the ith column of $Q$, to get an equivalent\n",
    "factorization with $R_{ii} > 0$. However this step is not needed in practice, since\n",
    "negative elements on the diagonal do not pose any problem in applications of the\n",
    "$QR$ factorization."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 67,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "4×4 Array{Float64,2}:\n",
       " -2.33103  -0.391209   1.48122   1.03352 \n",
       "  0.0      -1.93733    1.20418   0.192516\n",
       "  0.0       0.0       -2.02817   0.282745\n",
       "  0.0       0.0        0.0      -1.45234 "
      ]
     },
     "execution_count": 67,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A = randn(6,4);\n",
    "Q, R = qr(A);\n",
    "R"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 68,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "6×4 Array{Float64,2}:\n",
       " -0.524284    -0.393571    0.0235341   0.619977 \n",
       "  0.673821    -0.267557    0.12604    -0.0666858\n",
       " -0.420106     0.106823   -0.45678    -0.604267 \n",
       " -0.116714     0.133506    0.437192   -0.174386 \n",
       " -0.00913436  -0.857645   -0.12112    -0.330596 \n",
       "  0.284423     0.0934077  -0.75439     0.326097 "
      ]
     },
     "execution_count": 68,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Q = Matrix(Q)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 69,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1.2423769326346545e-15"
      ]
     },
     "execution_count": 69,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "norm(Q*R-A)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 70,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "4×4 Array{Float64,2}:\n",
       "  1.0           1.29462e-16  -4.52952e-17  -7.10347e-17\n",
       "  1.29462e-16   1.0          -4.09463e-18  -1.68385e-17\n",
       " -4.52952e-17  -4.09463e-18   1.0           2.43491e-16\n",
       " -7.10347e-17  -1.68385e-17   2.43491e-16   1.0        "
      ]
     },
     "execution_count": 70,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Q'*Q"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Julia 1.1.1",
   "language": "julia",
   "name": "julia-1.1"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "1.1.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
